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Abstract 
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reduced transition probabilities, phase-shifts, photo-disintegration cross sections, radiative capture 
cross sections, and astrophysical S-factors, for a two-body nuclear system. The code is based on a 
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reaction rates in numerous astrophysical scenarios. 
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PROGRAM SUMMARY 

1. Title of program: RADCAP (RADiative CApture) 

Computers: The code has been created on an IBM-PC, but also runs on UNIX ma- 
chines. 

Operating systems: WINDOWS or UNIX 
Program language used: Fortran- 77 

Memory required to execute with typical data: 8 Mbytes of RAM memory and 2 MB 
of hard disk space 

No. of bits in a word: 32 or 64 

Memory required for test run with typical data: 2 MB 

No. of lines in distributed program, including test data, etc.: 3054 

Distribution format: ASCII 

Keywords: Potential model; Photodissociation; Radiative capture; Astrophysical S- 
factors 

Nature of physical problem: The program calculates bound and continuum wavefunc- 
tions, phase-shifts and resonance widths, astrophysical S-factors, and other quantities 
of interest for direct capture reactions. 

Method of solution: Solves the radial Schrodinger equation for bound and for contin- 
uum states. First the eigenenergy is estimated by using the WKB method. Then, a 
Numerov integration is used outwardly and inwardly and a matching at the nuclear 
surface is done to obtain the energy and the bound state wavefunction with good accu- 
racy. The continuum states are obtained by a Runge-Kutta integration, matching the 
Coulomb wavefunctions at large distances outside the range of the nuclear potential. 

Typical running time: Almost all the CPU time is consumed by the solution of the 
radial Schrodinger equation. It is about 1 min on a 1GHz Intel P4-processor machine 
for a Woods-Saxon potential. 
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LONG WRITE-UP 



I. INTRODUCTION 



In astrophysically relevant nuclear reactions two opposite reaction mechanisms are of 
importance, compound-nucleus formation and direct reactions (for more details, see, e.g., 
Q]). At the low reaction energies occurring in primordial and stellar nucleosynthesis the 
direct mechanism cannot be neglected and can even be dominant. The reason for this 
behavior is that only a few levels exist for low excitations of the compound nucleus. 

In order to calculate the direct capture cross sections one needs to solve the many body 
problem for the bound and continuum states of relevance for the capture process. There 
are several levels of difficulty in attacking this problem. The simplest solution is based 
on a potential model to obtain single-particle energies and wavefunctions. In numerous 
situations this solution is good enough to obtain the cross sections within the accuracy 
required to reproduce the experiments. 

In this article a computer program is described which aims at calculating direct capture 
cross sections, based on a potential model. The program calculates bound and continuum 
wavefunctions, phase-shifts, energy location of resonances, as well as the particle-decay 
width, photodisintegration cross sections, radiative capture cross sections and astrophysical 
S-factors. The formalism for this model has been developed in Refs. 



II. BOUND STATES 



The computer code RADCAP calculates various quantities of interest for two-body fusion 
reactions of the type 

a + b — > c + 7, or a (6, 7)0. (1) 

The internal structure of the nuclei a and b is not taken into account. Thus, the states of the 
nucleus c is obtained by the solution of the Schrodinger equation for the relative motion of 
a and b in a nuclear + Coulomb potential. Particles a, b, and c have intrinsic spins labelled 
by I a, h and J, respectively. The corresponding magnetic substates are labelled by M a , Mj, 
and M. The orbital angular momentum for the relative motion of a + b is described by / and 
m. In most situations of interest, the particle b is a nucleon and a is a "core" nucleus. Thus 
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it is convenient to couple angular momenta as 1 + I&= j and j + I a = J, where J is called the 
channel spin. Below we also use the notation s, instead of I&, for the intrinsic spin of particle 
b. 

The bound state wavefunctions of c are specified by 

*JM (r) = ^y l JM , (2) 

where r is the relative coordinate of a and b, uf^ (r) is the radial wavefunction and y l JM is 
the spin-angle wavefunction 

y l jM= Yl (jmI a M a \JM) \jm) \I a M a ) , with \jm) = ^ Y lmi (?) X M b (3) 

m , M a mi , M b 

where XM b is the spinor wavefunction of particle b and {jmI a M a \JM) is a Clebsch-Gordan 
coefficient. 

The ground-state wavefunction is normalized so that 

CO 

J d 3 r |^ JM (r)| 2 = Jdr |«£ (r)| 2 = 1. (4) 
o 

The wavefunctions are calculated using a spin-orbit potential of the form 

V(r) = V (r) + V s (r) (l.s) + V c (r) (5) 

where V (r) and Vg{r) are the central and spin-orbit interaction, respectively, and Vc(r) is 
the Coulomb potential of a uniform distribution of charges: 

Vc( r ) — a b6 f° r r > Rc 

r 

Z a Z b e 2 / r 2 



where Z { is the charge number of nucleus i — a,b. 

One can use two kinds of approach to build up the potentials Vb(r) and Vs(r). In a 
Woods-Saxon parametrization they are given by 

h \ 2 1 d 



(7) 



K,(r) = V / (r), and V s (r) = - F 50 ( — ) - ^fs(r) 



with /j(r) 



1 + exp 



r-R, 
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The spin-orbit interaction in Eq. [7| is written in terms of the pion Compton wavelength, 
h/m n c = 1.414 fm. The parameters Vb, Vso, Ro, clq, Rsq, and 050 are adjusted so that the 
ground state energy E B (or the energy of an excited state) is reproduced. 

Alternatively, and perhaps more adequate for some situations, one can construct the 
potentials using a more microscopic approach. Among these models, the M3Y interaction is 
very popular. It has been shown to work quite reasonably for elastic and inelastic scattering 
of heavy ions at low and intermediate energy nuclear collisions 0, |(| . It has been applied to 
calculations of radiative capture cross sections with relative success (see, e.g., J?)). 

In its simplest form the M3Y interaction is given by two direct terms with different ranges, 
and an exchange term represented by a delta interaction: 

p-Pis p-fts 

t(s) = A— + B— + C5(s) , (8) 
Pis p 2 s 



where one of the possible set for these parameters is given by p, |6J A — 7999 MeV, B = 
-2134 MeV, C = -276 MeV fm 3 , (3 1 = A fm- 1 , and (3 2 = 2.5 fm' 1 . 

The central part of the potential is obtained by a folding of this interaction with the 
ground state densities, p a and p b , of the nuclei a and b: 

V M3Y (r) = X V M3Y (r) = X j d 3 n d 3 r 2 p a (n)p b (r 2 ) t(s) , (9) 

with s — |r + r 2 — ri|. Ao is a normalization factor which is close to unity. We assume that 
;he densities pi are spherically symmetric. The nuclear densities can be taken from e.g. Ref. 
I| for the charge matter densities. To obtain the matter density one can use the relation 

( r m) 1,/2 = \j ( r ch) — (0.85) 2 , where (r^) 1 ^ 2 and (r^) 1 ^ 2 are the charge and matter rms radii 

of the nucleus and the proton radius is taken as 0.85 fm. 

The spin-orbit part of the optical potential is parametrized as 

V s M3Y (r) = -A 50 (— V - ^V M3Y {r). (10) 
The bound-state wavefunctions are calculated by solving the radial Schrodinger equation 



h 



2 



2m ab 



df_ _ 1(1 + 1) 

dr 2 r 2 



u J u (r) + [V (r) + V c (r) + (s.l) V so (r)] (r) = E iU J tj (r) (11) 



where (s.l) = + 1) — /(/ + 1) — s(s + 1)] /2. This equation must satisfy the boundary 
conditions ufj (r = 0) = ufj (r = 00) = which is only possible for discrete energies E 
corresponding to the bound states of the nuclear + Coulomb potential. 



III. CONTINUUM STATES 



The continuum wavefunctions are calculated with the potential model as described above. 
The parameters are often not the same as the ones used for the bound states. The continuum 
states are now identified by the notation u El Ar), where the (continuous) energy E is related 
to the relative momentum k of the system a + b by E = h 2 k 2 /2m a b. 

The radial equation to be solved is the same as Eq. El but with the boundary conditions 
at infinity replaced by (see, e.g, Ref. 0) 



m ab 

UmAT > OO] 



if," (r) - S U H™ (r) e*™ (12) 



iai{E) 



where Su = exp [2iSu (E)], with 5u (E) being the nuclear phase-shift and ai(E) the 
Coulomb one, and 

fl I (±) (r) = G?,(r)±iF,(r) . (13) 

Fi and Gi are the regular and irregular Coulomb wavefunctions. If the particle b is not 
charged (e.g., a neutron) the Coulomb functions reduce to the usual spherical Bessel func- 
tions, ji (r) and rii (r). 

At a conveniently chosen large distance r = R, outside the range of the nuclear potential, 
one can define the logarithmic derivative 

( du J Elj /dr\ 

The phaseshifts 5u (E) are obtained by matching the logarithmic derivative with the asymp- 
totic value obtained with the Coulomb wavefunctions. This procedure yields 

_ G[- iF[ - a u (d - ig) 
lJ q + iFl-au (Gt + iFtY 1 ' 

where the primes mean derivation with respect to the radial coordinate at the position R. 

The continuum wavefunctions are normalized so as to satisfy the relation 

u J Elj \u J m/j ,) =5{E- E') 6jj,6 jr 6 w , (16) 
what means, in practice, that the continuum wavefunctions UEij{r) are normalized to 



— ^J2m a } 1 /^h 2 k e t5lJ sin(Ar + 5u) at large r. 

A resonance in a particular channel I J is characterized by 

dS u 



dE 2 



0, and 



r dE 



> . (17) 



j 



The single-particle width of the resonance can be calculated from 

/ \ _1 

2 / d5 u * 



1 IJ 



dE 



(18) 



E R 



IV. MULTIPOLE MATRIX ELEMENTS 

The operators for electric transitions of multipolarity \tc are given by (see, e.g. Ref. [ill ]) 

Oex^ = e x r A Y^ (?) , (19) 

where the effective charge, which takes into account the displacement of the center-of-mass, 
is 



e A = Z b e - 



m n 



m ( 



A 



For magnetic dipole transitions 

y 



O 




UN 



m n 



eM 



A 



m 2 a Z a m\Z h 



+ 



mi 



mi 



(20) 



(21) 



where Z M and are the spherical components of order \i (^t = —1, 0, 1) of the orbital and 
spin angular momentum (1 = — ir x V, and s = a/2) and gi are the gyromagnetic factors of 
particles a and b. The nuclear magneton is given by /xjv = ehjlm^c. 

The matrix element for the transition JoM — > JM, using the convention of Ref. 
is given by 



(JM J M > = (J M \fj,\JM) 



(J\\O EX \\J ) 



(22) 



V2J + 1 

From the single-particle wavefunctions one can calculate the reduced matrix elements 
(Ij \\Oex\\ hio) j- The subscript J is a reminder that the matrix element depends on the 
channel spin J, because one can use different potentials in the different channels. The 
reduced matrix element (J ||C.ea|| Jo) can be obtained from a standard formula of angular 
momentum algebra, e.g. Eq. (7.17) of Ref. One gets 



{J\\Oex\\Jo) 



-1 



>i + /a + J0 + A 



[(2J+1)(2J + 1)] 



1/2 



(Ij\\Oex\\IoJo)j- (23) 



j J I a 

Jo jo A 

To obtain (Ij \\Oex\\ hjo) j one needs the matrix element (ij ||r A Y\|| ^ojo)j for the spherical 



harmonics, e.g. Eq. (A2.23) of Ref. |l3j]. For l + I + A = even, the result is 

e A , -, ,/,,-!-/-;'!,-; /x ./n . ■ \ x , . i . / , \ ./ , , J 



(lj\\0 E x\\lojo)^ 



An 



J 



dr r A ujj (r) <° io (r) , (24) 



where we use here the notation k = \J1k + 1, and k = \Jk{k + 1). For l + 1 + A = odd, the 
reduced matrix element is null. 

Eqs. 1221 and 1221 can also be used for the magnetic dipole excitations. In comparison with 
the electric dipole transitions their cross sections are reduced by a factor of v 2 /c 2 , where v is 
the relative velocity of the a + b system. At very low energies, »<c, and the Ml transitions 
will be much smaller than the electric transitions. Only in the case of sharp resonances, the 
Ml transitions play a role, e.g. for the J = 1 + state in 8 B at = 630 keV above the 
proton separation threshold jj, [3| . However, the potential model apparently is not good 
in reproducing the Ml transition amplitudes |l5[. We only treat here the case in which the 



l jo) j and 



particle b is a nucleon. For that one needs the reduced matrix elements ( I j 

rrr 

(h ll^ll ^oio)j which are, e.g., given by Eqs. (A2.20) and (A2.19) of Ref. [13]. For I = l one 
obtains 



(Ij \\0 M i\\ Wo) t 



-n'-+'"+-'"+' J± jjA 3 J Ia 

47r 1 J Jo 1 



^o 



-J^ { l 0^o, '0+1/2 + (k + 1) Sj Q) l -l/2) + (-1)' 0+1/2 3 -^=S jQt l ±l/%8 jt k)T i/2 

l V2 



+ 9n 



1 

T 2 

'o 



(-1) 



Zo+l/2-jo 



~ x / -, \l +l/2-j JO r x 

JoOj, j - {-*-) 0j , l ±l/2<Jj, loTl/2 



+9a{-l 



)Ia+j0+J+1 ~ JIJa 




dr u J l3 (r) uf Q ° jo (r) 



(25) 



The spin g-factor is = 5.586 for the proton and = —3.826 for the neutron. The 
magnetic moment of the core nucleus is given by /i a = g a ^N- If I ^ la the magnetic dipole 
matrix element is zero. 
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V. THE ASTROPHYSICAL S-FACTOR 



The multipole strength, or response functions, for a particular partial wave, summed over 
final channel spins, is defined by 

dB (7rA; /qJq — -+ klj) _ |(fcJ||a A || Jo) | 2 
dk 2 J + 1 

= ^(2J+1) J 3 J h 1 \(klj\\O vX \\loj )j\ 2 , (26) 
j ( Jo jo A J 

where n = E, or M. 

If the matrix elements are independent of the channel spin, this sum reduces to the usual 
single-particle strength | (klj \\On\\\ Z jo) f / (2jo + !)• For transitions between the bound 



states the same formula as above can be used to obtain the reduced transition probability 

l klj ( r ) 

by the bound state wavefunction ujj 



by replacing the continuum wavefunctions ujL- (r) by the bound state wavefunction uL (r) 



That is, 

£?(7rA; ZtfoJo— tf.7) = (2J + l) J J J 4 1 \(lj \\O vX \\ l j )\ 2 . (27) 

[ Jo Jo A J 

For bound state to continuum transitions the total multipole strength is obtained by 
summing over all partial waves, 

dB (ttA) ^ dB (ttA; l j — ► klj) 



^ • (28) 



dE ^ dE 

ij 

The differential form of the response function in terms of the momentum E is a result of 
the normalization of the continuum waves according to Eq. 

The photo-absorption cross section for the reaction 7 + c — > a + c is given in terms of 
the response function by [l(3] 

< X) W - | 2 ^ (A / ,'; ( f ^ (29) 



A[(2A + 1)!!] 2 \h 2 k)\hc) dE 



where E 1 = E+\Eb\, with \Eb\ being the binding energy of the a + b system. For transitions 
between bound states, one has 



(30) 



where E{ (Ef) is the energy of the initial (final) state. 

The cross section for the radiative capture process a + b — > c + 7 can be obtained by 



detailed balance 



lOj . and one gets 



The total capture cross section cr nr is determined by the capture to all bound states with 
the single particle spectroscopic factors C 2 Si in the final nucleus 

a nT {E) = Y,{C 2 S) l a { Cl li {E). (32) 

Experimental information or detailed shell model calculations have to be performed to obtain 
the spectroscopic factors (C 2 S)i. For example, the code OXBASH [l6[ can used for this 
purpose. 

For charged particles the astrophysical S-factor for the direct capture from a continuum 
state to the bound state is defined as 

S {c) (E) = E a nT (E) exp [2^77 (E)} , with 77 (E) = Z a Z b e 2 /Hv, (33) 

where v is the relative velocity between a and b. 



VI. NUCLEAR REACTION RATES IN STELLAR ENVIRONMENTS 

The nuclear reaction rate, measuring the number of reactions per particle pair, a + b, 
per second in the stellar environment can be calculated from the nuclear cross section a 
for a given reaction by folding it with the velocity distribution of the particles involved. In 
most astrophysical applications the nuclei are in a thermalized plasma, yielding a Maxwell- 
Boltzmann velocity distribution. The astrophysical reaction rate R at a temperature T can 
then be written as 

R(T) = M , (34) 

1 + Oab 

where n« is the number density of the nuclear species i. The denominator takes care of the 
special case of two identical nuclei in the entrance channel. The quantity (crv) is given by 

1/2 



(av) 



(—) n^FWi r a ^ E exp (-ITr) dE ' ( 35 ) 
\7rm ab J {k B Ty/ 2 J \ k B T J 
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with &b the Boltzmann constant. 

The threshold behavior of radiative capture cross sections is fundamental in nuclear as- 
trophysics because of the small projectile energies in the thermonuclear region. For example, 
for neutron capture near the threshold the cross section can be written |lQj as 

7i —AkR Ima , or , 

a if = 7£ 1 [2 > I 36 ) 

k \a \ 

where «o is the logarithmic derivative for the s wave, given by Eq. Since «o is only 
weakly dependent on the projectile energy, one obtains for low energies the well-known 
1 / f-behavior. 

With increasing neutron energy higher partial waves with I > contribute more signifi- 
cantly to the radiative capture cross section. Thus the product av becomes a slowly varying 
function of the neutron velocity and one can expand this quantity in terms of v or \/E 
around zero energy: 

av = S {n) {0) + S {n) {0)VE + ^S {n) {0)E + ... . (37) 

The quantity S {n) (E) = av is the astrophysical S-factor for neutron-induced reactions and 
the dotted quantities represent derivatives with respect to E 1 ^ 2 , i.e., = 2\[E and 
g(v>) = 4E d ^ n) + 2 ^fj=r-. Notice that the above astrophysical S-factor for neutron-induced 
reactions is different from that for charged-particle induced reactions. In the astrophysi- 
cal S-factor for charged-particle induced reactions also the penetration factor through the 
Coulomb barrier has to be considered (Eq. l33|) . 

Inserting this into Eq. ESI we obtain for the reaction rate for neutron-induced reactions 

(av) = 5(0) + f^j 2 5(0)(A;bT)5 + ^S(0)k B T + .... (38) 

In most astrophysical neutron-induced reactions, neutron s-waves will dominate, result- 
ing in a cross section showing a 1/f-behavior (i.e., cr(E) oc 1/VE). In this case, the reaction 
rate will become independent of temperature, R = const. Therefore it will suffice to mea- 
sure the cross section at one temperature in order to calculate the rates for a wider range 
of temperatures. The rate can then be computed very easily by using 

R = (av) = (a) T VT = const , (39) 

with 



Vt 



2kT\ 



1/2 



m 



1 . (40) 
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The mean lifetime r n of a nucleus against neutron capture, i.e., the mean time between 
subsequent neutron captures is inversely proportional to the available number of neutrons 
n n and the reaction rate _R n7 : 

T » = • (41) 

If this time is shorter than the beta-decay half-life of the nucleus, it will be likely to capture 
a neutron before decaying (r-process). In this manner, more and more neutrons can be 
captured to build up nuclei along an isotopic chain until the beta-decay half-life of an 
isotope finally becomes shorter than r n . With the very high neutron densities encountered 
in several astrophysical scenarios, isotopes very far-off stability can be synthesized. 

For low | i?B|-values, e.g. for halo- nuclei, the simple 1/f-law does not apply anymore. A 
significant deviation can be observed if the neutron energy is of the order of the l-E^I-value. 
In this case the response function in Eq. EElcan be calculated analytically under simplifying 
assumptions (see Ref. For direct capture to weakly bound final states, the bound-state 
wave function Uij{r) decreases very slowly in the nuclear exterior, so that the contributions 
come predominantly from far outside the nuclear region, i.e., from the nuclear halo. For this 
asymptotic region the scattering and bound wave functions in Eq. 21 can be approximated 
by their asymptotic expressions neglecting the nuclear potential [jj 

ui(kr) oc ji(kr), u h (r) oc h^\ir]r) , 

where ji and h^' are the spherical Bessel, and the Hankel function of the first kind, re- 
spectively. The separation energy \E B \ in the exit channel is related to the parameter rj by 
\E B \ = h 2 r ] 2 /(2m ab ). 

Performing the calculations of the radial integrals in Eq. one readily obtains the 
energy dependence of the radiative capture cross section for halo nuclei For example, 

for a transition s — >p it becomes 



(rc) / x 1 (E+ 3|-Efi|) 2 , 

ff (W S ^ p)0C 7f E+\E B \ ' (42) 
while a transition p— >s has the energy dependence 



a 



(rc) 



If E \Eb\ the conventional energy dependence is recovered. From the above equations 
one obtains that the reaction rate is not constant (for s-wave capture) or proportional to T 
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(for p-wave capture) in the case of small | .Eel-values. These general analytical results can 
be used as a guide for interpreting the numerical calculations involving neutron halo nuclei. 

In the case of charged particles S(E) is expected to be a slowly varying function in energy 
for non-resonant nuclear reactions. In this case, S (E) can be expanded in a McLaurin series, 



S (E) = S (0) + S (0) E + -S (0) E 2 



(44) 



Using this expansion in Eq. ESI and approximating the product of the exponentials 
exp (— E/ksT) and exp \2tii] (E)} by a Gaussian centered at the energy E , Eq. ESI can 
be evaluated as jl| 



{av) 



m ab 



1/2 



(kT) 



S cS (E ) exp I — — 



3/2 



(45) 



with 



S eS (E )=S(0) 



5 S(Q) 

1 H 1 — 

12r S(0) 



En 



35E 
12r 



s(o) 

2S (0) 



El 



89E 2 
12r 



(46) 



The quantity E defines the effective mean energy for thermonuclear fusion reactions at 
a given temperature T, 



£ = 1-22 {Z 2 a Z 2 b m ah Tl) 1/2 keV , 
where Tq measures the temperature in 10 6 K. The quantities r and A are given by 



T 



3£ 



A 



(E kT) 



1/2 



(47) 



(48) 



kT' ~ VS 

An analytical insight of the cross sections and astrophysical S-factors for proton-halo 
nuclei can also be developed (see, e.g., Ref. jscj)- However, due to the Coulomb field 
the expressions become more complicated. The analytical formulas for direct capture cross 
sections involving neutron and proton halo nuclei are very useful to interpret the results 
obtained in a numerical calculation. 

For the case of resonances, where E T is the resonance energy, we can approximate a (E) 



10, 



2l|: 



by a Breit-Wigner resonance formula 

ixh 2 (2 J R + 1) 



°r{E) 



r r 



(49) 



2fiE (2J a + 1)(2J 6 + 1) (E r - Ef + (r tot /2) 2 ' 

where Jr, J a , and Jb are the spins of the resonance and the nuclei a and b, respectively, and 
the total width r to t is the sum of the particle decay partial width T p and the 7-ray partial 
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width T 7 . The particle partial width, or entrance channel width, T p can be expressed in 
terms of the sin gle- particle spectroscopic factor S and the single-particle width r s p of the 



resonance state 



22 



r p = C 2 S x r, p . , (50) 



where C is the isospin Clebsch-Gordan coefficient. The single-particle width r sp . can be 
calculated from the scattering phase shifts of a scattering potential with the potential pa- 
rameters being determined by matching the resonance energy (see Eq. I18jl . 

The gamma partial widths T 7 are calculated from the electromagnetic reduced transition 
probabilities B(Jj — > J/;L) which carry the nuclear structure information of the resonance 
states and the final bound states (^j). The reduced transition rates can be computed within 
the framework of the shell model. 

Most of the typical transitions are Ml or E2 transitions. For these the relations are 

r E2 [eV] = 8.13 x 1(T 7 [MeV] B(E2) [e 2 fm 4 ] (51) 

and 

r M i[eV] = 1.16 x 1(T 2 E% [MeV] B(M1) [fi 2 N ] . (52) 

For the case of narrow resonances, with width r -C E r , the Maxwellian exponent 
exp {—E/k-QT) can be taken out of the integral, and one finds 

3/2 



(<TV) 



fc^f) ^^) fl exp(-§), (53) 
where the resonance strength is defined by 



For broad resonances Eq. EH is usually calculated numerically. An interference term has 
to be added. The total capture cross section is then given by |24j 

a(E) = a m (E) + a T (E) + 2 [a ni (E)cr v (E)} 1/2 cos[5 R (E)} . (55) 

In this equation 5r(E) is the resonance phase shift. Close to a resonance, the phase shift 
approaches the value tt/2. Thus, close to a resonance one can use the expansion 

. , 7T , d8 

S R (E) ^--(E r -E)- 



R 
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Thus, using the definition given by Eq. one has 

r 

5 R (E) = arctan _ . (56) 

Only the contributions with the same angular momentum of the incoming wave interfere 
in Eq. ESI 



VII. COMPUTER PROGRAM AND USER'S MANUAL 

All nuclear quantities, either known from experiments or calculated from a model, as well 
as the conditions realized in the experiment, are explicitly specified as input parameters. 
The program RADCAP then computes the potentials, bound state energies, phase-shifts, 
transition probabilities, photo-dissociation cross sections and astrophysical .S-factors. 

The units used in the program are fm (femtometer) for distances and MeV for energies. 
The output cross sections are given in millibarns and the S-factors in eV.b. 

The program is very fast and does not require a complicated input. It asks the user the 
calculation one wants to perform. It is divided in 5 modules and one enters the following 
options when prompted on the screen: 

1 - for the calculation of M3Y potential, 

2 - for the calculation of energy and wavefunction of bound states, 

3 - for the calculation of reduced transition probabilities between bound states, 

4 - for the calculation of phase shifts and wavefunctions of continuum states, 

5 - for the calculation of astrophysical S-factors, response functions, photo-dissociation, 
and direct capture cross sections. 

For each option, a different subroutine is used: l=OMP_M3Y, 2=EIGEN, 3=BVALUE, 
4=CONT and 5=DICAP. 

The inputs can be commented by using the symbol "*" at the first position of an input 
line. 

Note that the angular momenta described in the text have the following correspondence 
in the program: I, l (L, LO); j, j (J, JO), I a , I b (s) (AIA, AIB), J, J (AICF, AIC). The 
program notation for the other variables are easy to recognize (see test cases below). 
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A. The M3Y potential 



To obtain the M3Y potential one selects the option 1. This calls the subroutine 
OMP_M3Y. The input file is named M3Y.INP. If the densities are not parametrized ei- 
ther by a Gaussian or by a Woods-Saxon function, one can enter them in the input file 
DENS.INP in rows of r x p a (r) x Pft(r). The first line of this input file should contain only 
the number of points in r. The function DENS_READ will read those densities and interpo- 
late them for use in the other routines. An example is calculation of the M3Y potential for 
the system p + 7 Be as follows. For the proton one can use a Gaussian density with radius 



parameter R = 0.7 fm. For 7 Be a Gaussian density parametrization can be used 8( with 
radius parameter R = 1.96. An appropriate input file is listed below. 

* ******** input of subroutine OMP.M3Y ******* 

* IOPT = Option for densities: = Gaussian or Woods-Saxon, 

* =1 densities entered in DENS.INP 

* NPNTS = number of points in the radial mesh (< 10000) 

* RMAX = maximum radius size (fm) ( < 250 fm) 

* IOPT NPTS RMAX 
100 10. 

* If IOPT = 0, enter density parameters: 

* Rl, Dl = Woods-Saxon form (radius and diffuseness) 

* R2, D2 = Same but for density of nucleus 2 

* For Gaussian densities, enter D1=0, or D2=0 

* Rl Dl R2 D2 

0. 7 0. 1.96 0. 

* Mass numbers 

* Al A2 

1. 7. 

The subroutine OMP_M3Y builds up the nuclear potential and calls the subroutine 
TWOFOLD which does all the work. It does the integration appearing in Eq. El The 
outputs will appear in files OMP.TXT and OMP.INP. The later one is for use as input of 
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FIG. 1: M3Y potential for the system p + Be. 

the M3Y potential by the other subroutines (if required). Figure 1 shows a plot of the 
potential obtained with this input. 

B. Eigenfunctions and energies 

The option 2 calls the subroutine EIGEN. If the real part of the potential is given as 
an input file OMP.INP (e.g. the one generated by the subroutine OMP_M3Y) it should be 
written in rows of r x V(r). The first line of this input file should contain only the number 
of points in r. The function OMP_READ will read and interpolate the potential for use 
in the other routines. The subroutine DERIVATIVE calculates its derivative to be used in 
the calculation of Eq. Let us assume we want to find the ground state of 8 B. The 2 + 
ground state of 8 B can be described as a p3/2 proton coupled to the 3/2~ ground state 7 Be. 
The subroutine POTENT builds up the potential. An example of the input file, named 
EIGEN. INP, which uses a Woods-Saxon potential, is shown as follows. 

* ******** input of subroutine EIGEN ******* 

* IOPT = option for potentials: 1 (2) for Woods-Saxon (M3Y) 

* NPNTS = no. of integration points in radial coordinate ( < 10000) 

* RMAX = maximum radius size ( < 250 fm) 
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* IOPT NPTS RMAX 
1 9999 250. 

* N_0 = nodes of the Wave Function (exclude origin) 

* JO = single-particle angular momentum 

* L0 = orbital angular momentum 

* N_0 JO L0 
1.5 1 

* If IOPT = 1, enter: 

* V0 = depth of central potential 

* VSO = depth of spin-orbit potential 

* R0 = radius of the potential 

* AA = diffuseness of the potential 

* RSO = radius of the spin-orbit potential 

* AAS = diffuseness of the spin-orbit potential 

* RC = Coulomb radius (usually, RC = R0) 

* 

* WS = V_0 f(r,R0,AA) - V_S0 (l.s) (r_0"2/r) d/dr f(r,RS0,AAS) 

* f(r,R0,a) = [ 1 + exp((r-R_0)/a) ]'(-l) 

* r_0 = 1.4138 fm is the Compton wavelength of the pion. 

* 

* V0 R0 AA VSO RSO AAS RC 
-44.658 2.391 0.52 -9.8 2.391 0.52 2.391 

* If IOPT = 2, or else (but not 1), enter FC, FSO and RC 

* (in this case, insert a '*' sign in above row, or delete it) 

* FC = multiplicative factor of central part of M3Y potential 

* FSO = multiplicative factor of spin-orbit part of M3Y potential 

* RC = Coulomb radius 

* FC FSO RC 

* 1.5 0.2 2.391 

* Zl, Z2 = charges of the nuclei 

* Al, A2 = masses of the nuclei (in nucleon mass units) 

* Zl Al Z2 A2 
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FIG. 2: The solid line shows the ground state wavefunction of 8 B in the potential model and the 
dotted line shows the real part of the wavefunction for the 1 + resonance in 8 B at 630 keV (see text 
for details). 

1. 1. 4. 7. 

In the example input shown above the potential parameters were chosen so as to reproduce 
the proton separation energy in 8 B which is equal to 0.136 MeV. If the M3Y potential was 
used, an input of the parameters Ao (FC), Xso (FSO), and Rc in Eqs. l6\ l9l and fTOl is needed. 
Note that this input line was commented, as we did not use it. 

The calculations are mainly done in the subroutine BOUNDWAVE which solves the 
Schrodinger equation for the bound-state problem. When Woods-Saxon potentials are used 
they are constructed in the routine POTENTIAL. 

The output of the wavefunction will be printed in EIGEN.TXT and GSWF.INP. The later 
is prepared for use as input wavefunction for the subroutine BVALUE (reduced transition 
probabilities), or the subroutine DICAP (direct capture subroutine). The solid line in Figure 
2 shows the ground state wavefunction of 8 B obtained with this input. 
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C. Reduced transition probabilities 

The option 3 calls the subroutine BVALUE which calculates reduced transition probabil- 
ities. To make and example with 8 B we artificially generate a p3/2, 1 + state, with excitation 
energy of 90 keV. This can be obtained by changing the WS potential input of EIGEN.INP 
to the values shown below 

* V0 R0 AA VSO RSO AAS RC 
-30.55 2.95 0.52 -8.53 2.95 0.52 2.95 

The output yields a state with energy of —0.05 MeV. 

The input file with the option 3 is to be given in BVALUE. INP. To calculate the reduced 
transition probability the input file should look like the example below. 

********* Input of subroutine BVALUE ******* 

* AIA = spin of the particle A (core) 

* AIB = intrinsic spin of the particle B 

* AIC = total angular momentum of the ground state of C = A + B 

* (channel spin) 

* JO = single particle angular momentum of B respective to A 

* L0 = relative orbital angular momentum of the ground state 

* AIA AIB AIC JO L0 
1.5 0.5 2 1..5 1 

* N_l = nodes of the excited state wave function (exclude origin) 

* J = single-particle angular momentum 

* L = orbital angular momentum 

* AICF = spin of the excited state after all angular momentum coupling 

* (channel spin) 

* N_l J L AICF 
1.5 1 1 

* JOPT = 1 (0) if final state angular momentum, AICF, is (is not) to be 

* summed over all possible values. If JOPT=l, AICF in the 

* previous line can be entered as any value. 

* JOPT 
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* Zl, Z2 = charges of the nuclei 

* Al, A2 = masses of the nuclei (in nucleon mass units) 

* Zl Al Z2 A2 
1. 1. 4. 7. 

* IOPT = option for potentials: 1 (2) for Woods-Saxon (M3Y) 

* Integration parameters for radial wavefunctions: 

* NPNTS = no. of integration points in radial coordinate ( < 10000) 

* RMAX = maximum radius size ( 250 fm) 

* IOPT NPTS RMAX 
1 9999 250. 

* V0 = depth of central potential 

* R0 = radius of the central potential 

* AA = diffuseness of the central potential 

* VS0 = depth of spin-orbit potential 

* RS0 = radius of the spin-orbit potential 

* AAS = diffuseness of the spin-orbit potential 

* RC = Coulomb radius (usually, RC = R0) 

* 

* WS = V_0 f(r,R0,AA) - V_S0 (l.s) (r_0"2/r) d/dr f(r,RS0,AAS) 

* f(r,R0,a) = [ 1 + exp((r-R_0)/a) ]"(-l) 

* r_0 = 1.4138 fm is the Compton wavelength of the pion. 

* 

* V0 R0 AA VS0 RS0 AAS RC 
-30.55 2.95 0.52 -8.53 2.95 0.52 2.95 

* If IOPT = 2, or else (but not 1), enter FC, FSO and RC: 

* (in this case, insert a '*' sign in above row, or delete it) 

* FC = multiplicative factor of central part of M3Y potential 

* FSO = multiplicative factor of spin-orbit part of M3Y potential 

* RC = Coulomb radius 

* FC FSO RC 

* 1.5 0.2 2.391 
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* MP = multipolarity: (Ml), 1 (El), 2 (E2) 

* MP 
2 

* GA = magnetic moment (in units of the nuclear magneton) of 

* particle A (core) 

* GB = magnetic moment of particle B 

* GA GB 
2.79 -1.7 

The output of this run yields B (E2; i — > /) = 3.76 e 2 fm 4 . The spectroscopic factors 
for the initial and final states are taken as the unity. If they are known one just multiply 
this result by their corresponding values. 

The bound state is calculated by the routine BOUNDWAVE and the 3-j and 6-j coeffi- 
cients are calculated in the routines THREEJ and SIXJ, respectively. 

D. Phase-shifts and resonances 

If one uses the option 4 the program will calculate the scattering phase-shifts for a given 
set of potential parameters and angular momentum quantum numbers for the continuum 
waves. For example, one might want to calculate the phase-shifts for the p+ 8 Be system in 
the energy interval E = - 3 MeV. The input file CONT.INP could be written as follows. 

* ******** Input of subroutine CONT ******* 

* IOPT = option for potentials: 1 (2) for Woods-Saxon (M3Y) 

* NPNTS = no. of integration points in radial coordinate ( < 10000) 

* RMAX = maximum radius size (< 250 fm) 

* NEPTS = number of points in energy ( < 1000) 

* IOPT NPNTS RMAX NEPTS 

1 9999 250. 200 

* V0 = depth of central potential 

* VS0 = depth of spin-orbit potential 

* R0 = radius of the potential 
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* AA = diffuseness of the potential 

* RSO = radius of the spin-orbit potential 

* AAS = diffuseness of the spin-orbit potential 

* RC = Coulomb radius (usually RC = RO) 

* 

* WS = V_0 f(r,R0,AA) - V_S0 (l.s) (r_(T2/r) d/dr f(r,RS0,AAS) 

* f(r,R0,A) = [ 1 + exp((r-R0)/a) ]"(-l) 

* r_0 = 1.4138 fm is the Compton wavelength of the pion. 

* 

* VO RO AA VSO RSO AAS RC 

-42.3 2.391 0.52 -9.8 2.391 0.52 2.391 

* If IOPT = 2, or else (but not 1), enter FC, FSO and RC: 

* (in this case, insert a '*' sign in above row, or delete it) 

* FC = multiplicative factor of central part of M3Y potential 

* FSO = multiplicative factor of spin-orbit part of M3Y potential 

* RC = Coulomb radius 

* FC FSO RC 

* 1.5 0.2 2.391 

* Zl, Z2 = charges of the nuclei 

* Al, A2 = masses of the nuclei (in nucleon mass units) 

* Zl Al Z2 A2 

1. 1. 4. 7. 

* EI, EF = initial energy, final energy 

* L, J = orbital angular momentum, angular momentum j (1+s) 

* EI EF L J 

0. 3. 1 1.5 

A run with this input file will show the presence of a sharp resonance at 631 keV with a 
width of approximately 50 keV. As for the case of bound states, the same resonance can be 
obtained with a different set of WS potential parameters, e.g. with the parameters shown 
below. 

* V0 R0 AA VSO RSO AAS RC 
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E [MeV] 



FIG. 3: Phase-shift (solid line) and its derivative (dashed line) for the p+ 7 Be system with the 
potential parameters described in the text. The 1 + resonance at 630 keV is observed. 



-28.65 2.95 0.52 -8.5 2.95 0.52 2.95 

The continuum states are calculated by the subroutine CONTWAVE and the Coulomb 
wavefunctions are calculated by the subroutine COULOMB. 

The phase-shifts and their derivatives with respect to energy are printed in the output 
file CONT.TXT. The Figure 3 shows these quantities for the test case above. 

The program also allows for the output of the continuum wavefunction for a given energy. 
The output of the wavefunction is printed in CWAVE.TXT. The real part of the continuum 
wave function of the 1 + resonance state at 630 keV is shown in Figure 2 (dotted line). 

E. Direct capture cross sections 

By choosing the option 5 the program calculates the direct capture cross sections and 
related quantities. The input file is DICAP.INP. The calculations are done in the subroutine 
DICAP. The output files are DICAP.TXT where the strength functions (in units of e 2 
fm 2A ), photodissociation cross sections (in mb), direct capture cross sections (in mb), and 
the astrophysical S-factors (in eV.b) are printed; DICAPL.TXT where the same output is 
printed, prepared for using in a plot program; and SFAC.TXT where the S-factor and its 
first and second derivatives with respect to the energy are printed. These can be used in 
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the calculation of the reaction rates by using Eqs. ^3 and EHJ 
An input example is presented below. 

* ******** Input of program DICAP ******* 

* IOPT = option for potentials: 1 (2) for Woods-Saxon (M3Y) 

* NPNTS = no. of integration points in radial coordinate ( < 10000) 

* RMAX = maximum radius size (< 250 fm). 

* NEPTS = number of points in energy ( < 1000) 

* IOPT NPNTS RMAX NEPTS 
1 9999 250. 200 

* N_0 = nodes of the ground state wave function 

* AIA = spin of the particle A (core) 

* AIB = intrinsic spin of the particle B 

* AIC = total angular momentum of the ground state of C = A + B 

* (channel spin) 

* JO = single-particle angular momentum 

* L0 = orbital angular momentum 

* EBOUND = binding energy of the ground state (absolute value) 

* N_0 AIA AIB AIC AJ0 L0 EBOUND 
1.5 0.5 2 1.5 1 0.14 

* JOPT = 1 (0) if final state angular momentum, AICF, is (is not) to be 

* summed over all possible values. If JOPT=l, AICF can be 

* entered as any value. 

* AICF = spin of the excited state after all angular momentum coupling 

* (channel spin) 

* JOPT AICF 

1 1. 

* Zl, Z2 = charges of the nuclei 

* Al, A2 = masses of the nuclei (in nucleon mass units) 

* Zl Al Z2 A2 
1. 1. 4. 7. 

* V0 = depth of central potential 
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* RO = radius of the central potential 

* AA = diffuseness of the central potential 

* VSO = depth of spin-orbit potential 

* RSO = radius of the spin-orbit potential 

* AAS = diffuseness of the spin-orbit potential 

* RC = Coulomb radius (usually, RC = RO) 

* 

* WS = V_0 f(r,R0,AA) - V_S0 (l.s) (r_(T2/r) d/dr f(r,RS0,AAS) 

* f(r,R0,a) = [ 1 + exp((r-R_0)/a) ]'(-l) 

* r_0 = 1.4138 fm is the Compton wavelength of the pion. 

* 

* VO RO AA VSO RSO AAS RC 

-44.658 2.391 0.52 -9.8 2.391 0.52 2.391 

* If IOPT = 2, or else (but not 1), enter FC, FSO and RC: 

* (in this case, insert a '*' sign in above row, or delete it) 

* FC = multiplicative factor of central part of M3Y potential 

* FSO = multiplicative factor of spin-orbit part of M3Y potential 

* RC = Coulomb radius 

* FC FSO RC 

* 1.5 0.2 2.391 

* EI,EF = initial relative energy, final relative energy 

* EI EF 
0. 3. 

* NS1,NP1,NP3,ND3,ND5,NF5,NF7 = (1) [0] for inclusion (no inclusion) 

* of sl/2, pl/2, p3/2, d3/2, d5/2, f5/2, and f7/2 partial waves 

* NS NP1 NP3 ND3 ND5 NF5 NF7 

10 110 

* MP = multipolarity: (Ml), 1 (El), 2 (E2) 

* SF = Spectroscopic factor 

* MP SF 

1 1. 

* GA = magnetic moment (in units of the nuclear magneton) of 
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* particle A (core) 

* GB = magnetic moment of particle B (proton, neutron, alpha, etc.) 

* GA GB 

-1.7 5.58 

Only the results for the astrophysical S-factor, S17, for the reaction p+ 8 B will be shown. 
They are plotted in Figure 4, together with the experimental data of several experiments 



25, 26, 27ll2a|. The first three set of data (MSU, GSI-1, and GSI-2) were obtained by using 

n n 

the Coulomb dissociation method |2£|. The other experimental results |2a| were obtained via 
a direct measurement. The dashed line shows the result of the calculated S-factor, obtained 
with the bound state wavefunction calculated with the same Woods-Saxon parameters as in 
the above input file. The dashed line represents the S-factor one obtains by changing the 
bound and continuum states using another set of Woods-Saxon potential parameters, which 
yields the same binding for 8 B, namely: 

* VO RO AA VSO RSO AAS RC 

-30.55 2.95 0.52 -8.53 2.95 0.52 2.95 



VIII. THINGS TO DO 

1 - Use the input files described above and reproduce the Figures 1-4. 

2 - Try to reproduce some of the radiative capture cross sections presented in the compi- 



lation of Ref. 
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3 - Show that for neutron halo nuclei the radiative capture cross sections follow the 
dependence described by equations H2 and 
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FIG. 4: Astrophysical S-factor for the reaction p+ 8 B. The data are the experimental points of 

nnnn 

recent experiments [25J, |26J, 1221 UM ■ The solid and dashed curves are results of calculations with 
different choices of the Woods-Saxon potential which reproduces the binding of 8 B. 
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